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SUMMARY 


Turbulence modeling for high-speed compressible flows is described and discussed. Starting with the 
compressible Navier-Stokes equations, methods of statistical averaging are described by means of which 
the Reynolds-averaged Navier-Stokes equations are developed. Unknown averages in these equations 
are approximated using various closure concepts. Zero-, one-, and two-equation eddy viscosity models, 
algebraic stress models, and Reynolds stress transport models are discussed. Computations of supersonic 
and hypersonic flows obtained using several of the models are discussed and compared with experimental 
results. Specific examples include attached boundary-layer flows, shock- wave boundary-layer interactions, 
and compressible shear layers. From these examples, conclusions regarding the status of modeling and 
recommendations for future studies are discussed. 


INTRODUCTION 


In this report we will discuss turbulence models that are used in numerical simulations of complex 
viscous flows. Although there are many applications and uses of turbulence models, we will restrict our 
attention primarily to high-speed compressible flows. The material covered constitutes a brief survey of the 
essential features of turbulence models and their status in applications, but does not include many details 
important in practice. For these, the reader is encouraged to consult the references. 

Turbulence models are necessary in numerical simulations because of the impracticality of computing 
all scales of turbulent motion. Since these scales compose a range many orders in magnitude, the computer 
storage required to resolve all scales is much larger than the storage capacity currently available on the most 
powerful computers. Even if computers did exist with the required capacity, the computational speed of 
current computers is too slow to handle all but the simplest of problems. Thus approximate methods, or 
models of turbulence, are introduced to simplify and make the computations practical. 

There are several approaches to turbulence modeling depending on how many of the turbulent scales 
are included in the modeling process. A more rigorous approach is to use subgrid-scale modeling (also 
known as large-eddy simulation) in which only turbulent eddies equal to or smaller than the numerical grid 
sizes are modeled. In this case the largest eddies are computed, and because they move and deform in time, 
the calculations are necessarily unsteady. This results in relatively large computing times and restricts the 
applicability of subgrid modeling to fundamental studies. 

A more practical approach is to model all the scales of turbulent motion. The equations solved in 
this case are the Reynolds-averaged Navier-Stokes equations and the numerical solutions, which represent 
long time averages of the flow variables, are usually steady in time. This is the approach described here. 

The report is organized into several main sections. It begins with a discussion of averaging proce- 
dures and the development of the Reynolds-averaged Navier-Stokes and related equations. After a brief 
discussion of the various types of turbulence models available, attention is directed to describing a rep- 
resentative sample of eddy viscosity models including explicit modifications to account for high speeds 


and compressibility. This is followed by a section on results in which representative computations are 
discussed and compared with experimental measurements. The paper concludes with a section on the cur- 
rent status of turbulence modeling for hypersonic flows with recommendations for future experiments and 
computation. 


5 

REYNOLDS- AVERAGED NAVIER-STOKES EQUATIONS 


The basic differential equations used in numerical simulations are the Reynolds-averaged Navier- 
Stokes equations. These equations are derived from the compressible Navier-Stokes equations by an aver- 
aging process that will be described shortly. The time-dependent, compressible Navier-Stokes equations 
are written as follows: 


Pt + (puj)j = 0 

( pudt + ( priiUj + = 0 * = 1 , 2,3 ( 1 ) 

(pE) t + (pEuj + u ivy + q } )j m 0 

where subscript notation has been used for partial derivatives, i.e., ( ) t = d/dt, ( ),, = d/dxi, and the 
summation convention is used for repeated indices. The molecular stress tensor and heat flux vector are 
expressed as 

(2) 

qj kTj 

In these equations, p is density, u< are Cartesian velocity components, E is total specific energy, T 
is temperature, h is enthalpy, e is internal energy, Pr is the Prandtl number, and the Stokes hypothesis is 
imposed. Assuming a perfect gas with constant specific heats, these variables are related as follows: 

P m (n- Upe e = C V T h = C P T 

7 = Cp/Cv E = e + u,-u,7 2 Pr = Cpp/n 


In most applications, the Sutherland relation is used for molecular viscosity, i.e., p = AT n /( B + T), 
where n, A, and B are constants that depend on the gas. 


To derive the Reynolds-averaged Navier-Stokes equations an averaging operation is defined as 
follows: 

1 r t+T 

= 2 fJ,-T piXi ' S)dS 


p(Xj,t) 


(4) 


where 2T is the averaging interval, which is assumed to be large compared with the energy containing 
turbulent time scales, but small compared with the time scale of the mean or average motion. The mean 
density, p, in this sense is a slowly varying function of time. Although p depends on the averaging interval, 
it is tacitly assumed that a range of values for T exists for which p is practically independent of T and it is 
this range that is applicable in the averaging operation. 
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An alternate form of averaging which may be used in place of time-averaging is ensemble-averaging in 
which the averaging is performed over a large number of records or experiments. In this case the difficulty 
of selecting a time-averaging interval is not present and therefore this form of averaging is superior in many 
respects to time-averaging. Because results obtained using the two averaging forms are identical, we will 
continue to use the term time-averaging even though the preferred form is ensemble-averaging. 

The fluctuating density, p', is defined as the difference between the density and its average value, i.e., 

p' = p-p (5) 

Averages and fluctuating quantities for other variables such as p and u, are defined similarly. Although 
not strictly true unless T — ► oo, we assume here that = 0 and = p, which is consistent with ensemble- 
averaging. 

The time-averaged Navier-Stokes equations are obtained by averaging equation (1). The result in- 
volves averages such as pu,u ; and p/iu ; , which can be split into averages of mean and fluctuation quanti- 
ties, e.g., 

puiuj = (p + p')(u< + u’i)(uj + u'j ) 

= p UiUj + pu'iu'j + Uip'u'j + Ujp'u'i + p/u'iu'j (6) 

For incompressible flows, where p' = 0 , the last three terms in the above equation are absent. From 
this it is apparent that the compressible averaged equations will contain many more terms than the incom- 
pressible averaged equations. For this reason, an alternative form of averaging for velocity and energy 
variables has been developed, which leads to a form of the compressible averaged equations that is almost 
identical to the incompressible form. This is called mass-averaging (or Favre-averaging) in which the mean 
and fluctuating velocities and enthalpies are defined as follows: 

u, = pu7/p, h = ph/p 

u'i — Ui — u» , h" = h-h {) 

It is important to note that averages of fluctuating quantities are no longer zero, but finite, i.e., 

< = -pX/p , F = -W/p (8) 

but that mass-weighted averages of u" and h" are zero, i.e., 

P< = (p+p')< = 0 , pF = (p+p')F = 0 (9) 

By introducing mean and fluctuating quantities into equation (1) and averaging, we obtain the mass- 
averaged form of the compressible Navier-Stokes equations. These equations are written below. For sim- 
plicity, the bar and tilde notations have been omitted from averaged variables. 


Pt + ( puj)j = 0 

(puj)t + (pu,u ; + (Jij)j = 0 i = 1,2,3 (10) 

(pE) t + ( pEuj + UjCij + qj)j = 0 

/T — i , T 

®\j + Qj Qj Qj 
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where off , of } , etc., are defined as 


Reynolds stress tensor: 

G ij 

= pu'lu'l = (p + p'Wlu'l 

Molecular stress tensor: 


= p8{j /i(uij + Uj t i * 

Reynolds heat-flux vector: 


= ph"u'j = (p + p')h"u"j 

Molecular heat-flux vector: 

«r 

= -g-hj 
Pr J 

Total energy: 

B 

= e + k + UkUk/2 

Turbulent kinetic energy: 

k 

= pUkUk/2p 


( 11 ) 


In equation (11) it is assumed that p is independent of time in the averages leading to off, and qf . 

The goal of turbulence modeling is to relate the Reynolds stresses and heat fluxes to known mean- 
flow quantities such as velocity and temperature. This can be done in various ways which lead to different 
types of turbulence models. If the Reynolds stresses and heat fluxes are related algebraically to the mean- 
flow variables, the corresponding models are called algebraic stress models. The most important and 
simplest subclass of these models are eddy viscosity models which relate Reynolds stresses to strain rates 
(or velocity derivatives) in a manner identical with molecular stresses. Eddy viscosity models will be the 
primary focus of this paper. 

The simplest eddy viscosity models are the zero-equation models in which the eddy viscosity is mod- 
eled algebraically in terms of flow geometry and mean flow variables. More complicated turbulence mod- 
els have been developed in which the Reynolds stresses are defined by field equations.These equations 
are derived by manipulating the Navier-Stokes equations for mean and fluctuating quantities (ref. 1). The 
resulting equation for the Reynolds stress tensor is given below 


Reynolds stress equation: 

<°£)t 

Production tensor: 

pPik 

Dissipation tensor: 

ptik 

Reynolds flux tensor: 

Qikj 

Fluctuating stress: 

„F 

a H 


+ ( of k Uj + Qikj) j = p(Pik - c,*) 
= -(oyUkj + OkjUij) 

= -{ofjUkj + ajttij) 

= (puX'u" + a(}K + offjU'l) 

2 

— SijP /i( U i,j + Uj,i ~6»yUt ( fc) 


( 12 ) 


In these equations the fluctuating stress tensor, ofi, is interpreted to include both mean and fluctuating 
quantities. The dissipation tensor, pe lfc , contains both mean and fluctuating pressures and velocities and 
in this sense is a more general or extended definition than the conventional ones, which contain only 
fluctuating velocities. The production tensor is given directly in terms of the Reynolds stresses and mean 
velocity components, and thus requires no modeling. However, the dissipation and Reynolds flux tensors 
involve unknown averages and must be modeled. 
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A simplified form of the Reynolds stress equation is obtained by taking its trace and is called the 
turbulent kinetic energy equation, or TKE equation. This equation is written below 




TKE equation: 

Production: 

Dissipation: 

TKE flux: 


(pk) t + ( pkuj + q kj )j = p(P - e) 
pP = -aJjUij 

pe = -afju'lj 
Qkj = jP«u'j + afjU'l 


(13) 


The TKE equation forms the basis of several classes of turbulence models, including the one- and 
two-equation models, and the algebraic stress models. In these models, the (square of the) velocity scale 
of turbulence is given by the TKE and a length scale of turbulence is given either algebraically or in terms 
of a field variable governed by an equation similar to the TKE equation. 


EDDY VISCOSITY MODELS 


Eddy viscosity models are the simplest turbulence models in the sense that they model turbulent 
stresses and fluxes by analogy to molecular stresses and fluxes. This approach is generally referred to as 
the Boussinesq approximation. The models may be expressed in terms of an eddy viscosity function, fi T , 
and a turbulent Prandtl number, Pr T , as follows: 


of/ - pu'iu'j = jfik - p T ( U ij + «/,( - tu.*) 




Pr T =Sw 

K'p 


(14) 


where k t is the turbulent conductivity. With these models, the whole problem of modeling is reduced to 
defining the eddy viscosity and turbulent Prandtl number. The turbulent Prandtl number is usually assumed 
to be a constant of the order of unity, but it may vary between classes of problems. (It is normally set equal 
to 0.9 for boundary-layer problems.) The eddy viscosity function may be expressed in terms of length and 
velocity scale functions, l and q, as follows: 


pt = pk 


(15) 


The way I and q are determined defines the type of eddy viscosity model to be used. If l and q are 
determined algebraically from mean flow data, the models are referred to as zero-equation models. If l is 
determined algebraically, but q is determined from a field equation such as the TKE equation, i.e., equation 
(13), the model is referred to as a one-equation model. If both l and q are determined from field equations, 
the resulting model is called a two-equation model. For this report we will discuss zero- and two-equation 
models. 


5 



Eddy viscosity models were developed originally for incompressible flows and only later were ex- 
tended to compressible flows. Aside from the use of mass-averaging instead of time-averaging, there is 
very little difference in form between the two types of models. This is because in the initial investigations, 
which were restricted to attached transonic and moderately supersonic flows, it was found that the incom- 
pressible forms were quite satisfactory. As we shall see, the extensions to higher-speed flows in some 
simple cases are satisfactory, but other more complex cases require specific corrections for compressibility 
effects. 


Zero-Equation Models 


Zero-equation models are the simplest of eddy viscosity models in the sense that they do not make use 
of additional field equations. In this section we will discuss two widely used models that are representative 
of most other zero-equation models and that will be discussed in the results section. Unless otherwise 
stated, it is assumed that these models are applied at solid walls using no-slip boundary conditions 

Cebeci-Smith Model. - The model described here is a simplified version of the model described in 
reference 2. It is a two-layer model that uses Prandtl’s mixing length model (ref. 3) for the inner layer and 
Clauser’s model (ref. 4) for the outer layer. The model is expressed as follows: 

Mr - mtn(/irj,/iro) 

Mtv = pl 2 s , pto - 0 0168 p8*u e /I 

l = 0.4 yd, d = 1 — exp(— y+/A+) (16) 

Ut = 'Jtw/Pu , y + = u T y/u v 
6* = f\ 1 - u/u e )dy , 1=1 + 5 .5 (y/8) 6 

the strain-rate parameter, s, is usually taken to be the shearing strain , |u v + v x \, for two-dimensional 
problems. The operation of taking the minimum in equation (16) is interpreted to mean using the inner 
eddy viscosity, Mri. until it first becomes larger than the outer eddy viscosity, pro> beyond which point 
the outer formula is used exclusively. Parameter I is Klebanoff’s intermittency factor, u T , is the friction 
velocity, r w , is the wall shear, u w = p w /p w , and the subscript w indicates wall values. The nondimensional 
parameter A* , from van Driest (ref. 5) generally depends on the streamwise pressure gradient and surface 
blowing and roughness characteristics (ref. 2). For boundary layers with smooth solid walls and zero or 
small pressure gradients, A + is a constant, i.e., 


A + = 26 


A more complicated and general version of this model including transition modeling terms is given 
in reference 2. 

Baldwin-Lomax Model. - The Baldwin-Lomax (B-L) turbulence model (ref. 6) is similar to the 
Cebeci-Smith (C-S) model, but it incorporates features that make it more advantageous for complex 
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two- and three-dimensional flows. It is similar to the (C-S) model in that it uses nearly the same inner 
model, but it differs with respect to the outer model. This model can be expressed as follows: 

= Tnin(/ir/,/iro) 

Mn = pi 2 s , s = y/ujiuTi 

/J>TO = 0 .027 Pilmax nXJTlf F max > 0 .25 Ujj/ F max ) / 1 (17) 

F = ysd, Ud ~ | U | max |u) TO ,„ , 8 = l/moi/0 -3 

where y max is the outermost value of y in the boundary layer where F has a local maximum, F max . In these 
formulas l is the Prandtl mixing length given by equation (16), d is the Van Driest damping factor, and I 
is the intermittency factor of equation (16) in which 8 is replaced by y max / 0 .3 as indicated. 

It should be noted that with the B-L model, in contrast with the C-S model, the strain-rate function, 
s, is defined as the magnitude of the vorticity vector and not the shearing strain. This makes the model 
directly applicable to three-dimensional problems where an invariant shearing strain is not well defined. 

A basic advantage of this model over the C-S model is a result of how the outer model is defined. 
Referring to equation (16) we see that the C-S model requires both the displacement thickness, 8*, and the 
boundary-layer thickness, 8, which is used in the intermittency function. Both of these thickness parameters 
frequently are not well defined and are difficult to compute, especially for separated flows. The advantage 
of the B-L model is that it uses a length scale, y^s, which is well defined and easily computed for a wide 
class of flows. This does not necessarily mean that the B-L model is superior to the C-S model on a physical 
basis, but it does mean that it is more convenient on a numerical basis. 

From a physical standpoint, it has been found that the C-S and B-L models give similar predictions 
of both attached and separated boundary-layer flows for low to moderate supersonic flows. Predictions 
of attached flows are usually in good agreement with experiments, but predictions of separated flows are 
frequently deficient. At hypersonic speeds, the models also tend to give similar predictions, although there 
is some evidence that the B-L model may be more sensitive to Mach number than the C-S model is. 

The procedure of applying no-slip boundary conditions is frequently referred to as the integration-to- 
the-wall procedure. To be applicable, the numerical mesh spacing normal to the wall must be chosen such 
that the value of y + at the first point off the wall is of the order of unity, placing it well within the viscous 
sublayer. 

For some numerical algorithms, such as explicit methods, the procedure of integrating to the wall 
has detrimental effects on numerical stability because of the fine mesh spacing required. In this case an 
alternate approach called the law-of-the-wall procedure, or wall-function method, is used. In this approach 
a slip-type boundary condition, based on the logarithmic velocity law of turbulent boundary layers, is used. 
This law can be written (for incompressible flow) as 

u = — In Ey+ 

K 

k = 0.4, £7 = 9 .128 (18) 
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To apply this method, within the context of time-marching Navier-Stokes solvers, the above formula 
for velocity is solved (by a Newton-Raphson procedure) for the friction velocity, u T , using for u and y their 
values at the first point off the wall. The values of u w and p w are obtained by extrapolating the temperature 
to the wall. Once u r is determined, the wall shear stress, r w , is obtained, which is then used directly in the 
boundary condition for wall velocity. More complicated formulas have been developed for compressible 
flows and flows involving wall heat transfer. 7 


Two-Equation Models 

Zero-equation models are well adapted to simple attached flows where a single well defined shear 
layer is easily identified. There are many complex flows where this is not the case, however, and use 
of zero-equation models becomes difficult or unwieldy. Examples of such flows include separated flows 
behind bluff bodies and multiple intersecting shear layers. In these cases it is difficult to define appropriate 
velocity and length scales because several such scales are usually present in the flow. For this reason, more 
advanced models have been developed in which the velocity and length scales are determined from field 
equations. These are the two-equation eddy viscosity models. 

The prototype field equation for the two-equation models is the turbulent kinetic energy equation, 
equation (13). In order to use this equation, averages or correlations in the dissipation and TKE flux terms 
must be modeled in terms of known or mean-flow quantities. It is beyond the scope of this paper to explain 
in detail how these terms are modeled. Instead, we will simply discuss the results of the modeling. 

There are essentially two terms in the TKE equation that must be modeled. These are the TKE flux, 
q ki , and the dissipation, c. For the TKE flux, a gradient-diffusion approximation is used, i.e., 



where Pr k is a modeling constant (Prandtl number) of the order of unity. 

The absolute dissipation rate, e, is obtained from a separate field equation similar to the TKE equation 
given below. The velocity and length scales, and the eddy viscosity, are expressed in terms of k and e as 
follows 


Pt = C^fpql = Cufpk/u = C^fpk 2 /^ ( 20 ) 

q = Vk , l = q/u , u> - e/k 

where / is a damping function analogous to the van Driest damper, equation (16), is a modeling con- 
stant, and w is the specific dissipation rate. With these approximations, the TKE equation can be written 

( pk) t + ( pkuj + q k j)j ~ p( P - c) (21) 

_ 2 

pP = -aJjUij = prS - jpkD ( 22 ) 

2 

S = ( Uij + Uj ti )uij - ju 2 k k , D = 
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In these equations, P is the turbulent production that is reexpressed in terms of the eddy viscosity and the 
strain invariants S and D. For incompressible flows the dilatation, D = u k<k , is zero, but for compressible 
flows this term is nonzero and can be an important modeling term in some cases. 

The term on the right-hand side of equation (21) is a turbulence source function that may be expressed 
in terms of a nondimensional source function, /ijt, using equations (20) and (22), i.e.. 


p(F-e) = h k puk 

h -r f S 2D 1 
hk ~ ~ ~ 1 


ur 


3 u 


(23) 


Equation (21) is the prototype field equation for all two-equation models. With general two-equation 
models, the variables k and e (or w) of equation (20) are expressed in terms of two auxiliary variables, 
si , and 52, each of which is governed by field equations similar to equation (21). This general form of a 
two-equation model can be written as follows: 

k = k(si,s 2 ), € = e(si,s 2 ) 

( pSi)t + (pSiUj + Qij)j = hipwsi *- 1,2 

hi = Ci\ ( C „/4 - T~) ~ C a ( 24 ) 

\ w z 3 w/ 



In these equations there is no summation on the index i. 

The eddy viscosity damping function, /, is usually expressed in terms of a turbulence Reynolds num- 
ber, R t , which in turn is written in terms of k and e (or u>). Typical expressions for / and Rr are 

/ = 1 - exp( - cxR t ) , R t = — (25) 

pu> 

where a is a constant. 

For fully developed turbulent flows, the turbulent Reynolds number becomes very large and the damp- 
ing function approaches unity. On the other hand, at low Reynolds numbers (e.g., in laminar regions or the 
viscous sublayer) / goes to zero. In general, the variables Cn and C, 2 of equation (24) are also functions 
of the turbulence Reynolds number, analogous to equation (25), although in some cases they may involve 
additional terms. The Prandtl numbers Pr \ , Pr 2 are usually taken to be constant. At large values of the 
turbulence Reynolds number the variables Cn and C, 2 generally approach constant values along with the 
damper /. 
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1. k - e Model 


One of the most widely used two-equation models is the k — e model originated by Launder and 
Spalding (ref. 8). In this model si = k and 32 = e. The high Reynolds number form of the constants Cn 
and Ca, as well as the other constants in equation (24), are given by 

<^ = 0.09, Cn = Cn = f = Pr\ = 1 

C 21 = 1.45, 022 = 1.92, Pr 2 = 1.3 (26) 

These constants have been obtained by comparing solutions of the governing equations with experi- 
mental results. For example, the constants Cn and C 12 come directly from the TKE equation. The constant 
C 22 is determined from experiments on the decay of isotropic turbulence in which case all production and 
diffusion terms are absent from the equations, and an exact solution is easily obtained. The other constants 
are determined by obtaining approximate solutions for the wall region of equilibrium boundary layers, 
where P = c, and the logarithmic law, equation (18), is applicable, and by numerically optimizing free 
shear flow solutions. 

The high Reynolds number form of the k — e model described above is applicable to fully developed 
turbulent flows and does not apply to the viscous sublayer. For such applications, the molecular viscosity 
is much smaller than the turbulent viscosity and usually is neglected in the diffusion fluxes. In these cases, 
however, special slip-type boundary conditions based on equation (18) must be applied to the velocity and 
turbulence variables because the first numerical grid point must be taken well outside the viscous sublayer 
*(in the fully turbulent region) and no-slip conditions are inappropriate. This approach has been followed by 
Launder and Spalding (ref. 8) and others. The generalization to compressible flow is described by Viegas, 
Rubesin and Horstman (ref. 7). Although this approach is convenient for many problems it is not easily 
adapted to low Reynolds number flows where transition phenomena are important. In these cases, a more 
general low Reynolds number form of the model must be used in which Cn, Cn, etc., depend on R T . 
Several such models have been developed, including those by Jones and Launder (ref. 9), Chien (ref. 10), 
and Wilcox and Rubesin (ref. 11). Because the formulas defining these models are relatively complicated, 
we will not give them here. Instead, we will describe an alternative low Reynolds number, two-equation 
model that is given below. 

2a. q - w Model a 

The q — w model was developed to overcome numerical stability problems encountered with several 
low Reynolds number, two-equation models (refs. 9, 10, and 11). A discussion of these problems, and the 
development of the q — w model, is given in references 12, 13, and 14. For this model, the variables si 
and 3 2 of equation (24) are taken as 



s 2 = u - e/k 


(27) 


The parameters and constants in the equations are given by the following relations: 
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= 0 .09 , / = 1 — exp( —0 .02 qy/v) , Pn = Pr 2 = 2 
Cn = Cn = 0 .5 , C 2 i=0.055 + 0.5/ , C 22 = 0.833 , C 23 = | (28) 

h 2 = C 2 1 — Cm — 'j — C 22 

\ ur w/ 

Numerical boundary conditions to be applied with this model at solid walls are given by 
u = v = q = u) y = 0. 

2b. q — lo Model b: Compressibility correction 

The previous model was tested on an oblique shock-wave boundary-layer interaction flow for a sep- 
arated case, but it failed to predict any separation (ref. 15). In reference 15, a correction to the model 
was introduced that led to substantially improved predictions. This modification was based on the work 
of Morel and Monsour (ref. 16) who observed that in a uniaxial compression, the standard k — e model 
predicts that the turbulence length scale should increase, which runs counter to the physical expectation 
that it should decrease. Arguing that the product of the density and the turbulent length scale should remain 
constant in a uniaxial compression, they derived a correction to the source term of the e equation. Trans- 
lated to the w equation, this modification results in a new value of the constant multiplying the dilatation 
term, i.e., 

Cn = 2.4 (29) 

In reference 17, Vandromme proposed a compressibility modification with some similarities to the 
modification described here. His modification was based on earlier work by Rubesin and included density 
gradient terms as well as dilatation terms. Results using this model will be reported in the section on 
compressible shear layers. 

2c. q — w Model c: Heat transfer correction 

The previously described correction for compressibility improves the predictions of pressure distri- 
bution and separation, but the surface heat transfer remains relatively unaffected and too high in the region 
of reattachment. To remedy this difficulty, a modification or constraint on the turbulent length scale was 
imposed, following the work of Monsour reported in Kline, Cantwell, and Lilley (ref. 18). In this correc- 
tion, an upper bound is placed on the length scale appearing in the eddy viscosity such that it can never be 
greater than a constant times the Prandtl length scale in the wall region. The result is 

£= min(2Ay, q/w) (30) 


This correction generally does not change the predictions of the turbulence model exept near a reat- 
tachment point, and, to a lesser extent, near a separation point. This occurs because in equilibrium or 
attached flows, the turbulent length scale q/w is approximately equal to 2.4y in the wall region. This 
model will be referred to as the q — u model c. 
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RESULTS 


The status of modeling for high-speed flows will now be described by comparing predictions with the 
results of experiments. Experimentation plays an important role in the development of turbulence models 
by providing data on the flow physics required to substantiate modeling assumptions and in verifying the 
performance of models in testing. Wherever possible, experimental data from low- and high-speed flows 
will be contrasted to illustrate similarities and differences. Emphasis will be on attached flows, shock- wave 
boundary-layer interaction flows, and shear layers. References 19, 20, and 21 cite data that have been used 
to evaluate turbulence models for aerodynamic flow predictions. 


Attached Boundary-Layer Flows 

Modeling for hypersonic attached flows is more mature than for the other flows we consider. Eddy 
viscosity models perform reasonably well, as our examples will show. This fact may not be surprising 
because the modeling has been founded on a rather substantial experimental data base used together with 
knowledge regarding the behavior of incompressible flows. 

Figure 1 shows a composite sketch of a turbulent boundary layer constructed from a substantial incom- 
pressible data base. Velocity profile data can be collapsed onto a single curve using the friction velocity, 
u T , as a scaling parameter. Regions of the viscous sublayer, the logarithmic region, and the outer layer are 
depicted. The viscous sublayer is the region where molecular viscosity is important. It consists of a lami- 
nar sublayer region and a buffer region that blends with the logarithmic turbulent region. The logarithmic 
region is characteristic of all turbulent boundary layers and can be expressed as a function of the Reynolds 
number based on the friction velocity, or y + . The outer region, which actually begins quite close to the 
wall (y/6 between 0.1 and 0.2), is characterized by a wake-like region whose shape and thickness depend 
on the pressure gradient imposed by the outer inviscid flow field and the Reynolds number. 

At the high speeds associated with supersonic and hypersonic Mach numbers, similar experimental 
observations have been made. In these cases, however, it is necessary to introduce a compressiblity trans- 
formation (ref. 22) to adjust the profiles appropriately. Figure 2 shows the transformed velocity profiles 
taken in a very high Mach number helium wind tunnel. It should be noted that at very high Mach numbers, 
the pressure gradient must balance turbulent normal stresses arising form the normal momentum balance. 
Also, the sublayer becomes thicker as the Mach number increases. 

Representing turbulent velocity profiles using log-law variables enabled integration of the mean mo- 
mentum equation to determine such quantities as skin friction and heat transfer. But, with the advent of 
finite difference methods for solving the boundary-layer equations, the development of mixing length and 
eddy viscosity models was facilitated when it was experimentally observed that the shear stress distribution 
across a boundary layer changed little because of compressibility. Figure 3, taken from Sandbom (ref. 23), 
shows compressible data up to Mach 7, compared with similar data representative of incompressible flows. 
Earlier, Maise and McDonald (ref. 24), using a similar approach with adiabatic wall temperature data, 
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showed that the mixing length and eddy viscosity, scaled by the boundary-layer thickness and the incom- 
pressible displacement thickness, respectively, were essentially independent of compressibility effects up 
to a Mach number of 5. See figures 4a and 4b. 

The effects of compressibility and wall temperature on skin friction are shown in figures 5a and 5b. 
The solid line is the van Driest correlation based on the Karman-Schoener incompressible friction law and 
represents the available skin friction data to within 10% for the adiabatic wall data and to within 20% for 
the data with heat transfer. See Hopkins and Inouye (ref. 25). Computations using the boundary-layer 
equations are compared with the data in the figure. Aside from showing that eddy viscosity models predict 
the correct influence of compressibility on skin friction (fig. 5a), several other conclusions can be reached. 
The choice of mass-averaging or time-averaging has no significant effect on the predicted results. The zero- 
equation C-S model reproduces the van Driest result somewhat more accurately than do the other models; 
thus this model would have to be the choice for prediction, considering its simplicities. The effects of heat 
transfer are illustrated in figure 5b where a two-equation model prediction is compared with the van Driest 
correlation for M=5. The result, which is typical of most eddy viscosity predictions, deviates from the van 
Driest variation as total temperature ratio decreases and points to a caution regarding accurate prediction 
of cool-wall heat transfer trends, although the data are considerably scattered in these cases. 

Shang (ref. 26) extended computations using the C-S model to higher Mach numbers. He incorporated 
the normal momentum equation to account for nonzero normal pressure gradients and, more importantly, 
accounted for triple correlations involving density fluctuations usually omitted at lower Mach numbers. 
Results are shown in figures 6 and 7. Data and computations from two models, one with density fluctuation 
terms and one without density fluctuation terms, are compared. The inclusion of these terms affects the 
heat-transfer predictions somewhat more than the skin friction, but either approach produces reasonably 
accurate results, considering the uncertainties in the data. 

It is interesting to note that at lower Reynolds numbers the data tend to be underpredicted, particularly 
the heat transfer. Such results are common because boundary layer transition influences the region encom- 
passed by the low Reynolds number. These influences also tend to affect data correlation and may explain 
why there is more scatter in the cold-wall data around the van Driest predictions at higher Mach numbers 
where transition lengths are substantial. Figure 7, from Shang (ref. 26), shows skin-friction measurements 
for low and high Reynolds numbers compared with computations obtained with and without accounting 
for density fluctuation effects. At high Reynolds numbers, where the turbulent flow is fully developed, 
the computations compare reasonably well with the data, although it is difficult to conclude whether it is 
necessary to include the density fluctuation effects because of the data scatter. At low Reynolds numbers, 
the data are scattered for all Mach numbers, but this is especially pronounced at the highest Mach num- 
bers, probably because of transition effects. The computations show poorest agreement at the high Mach 
numbers, so a cautionary note is made for this regime. 

The effects of low Reynolds numbers can be accounted for at low Mach numbers approximately by 
modifying either the maximum mixing length or the outer eddy viscosity used in the model formulations. 
(See McDonald, ref. 27). Bushnell (refs. 28, 29, and 30) investigated the low Reynolds number problem 
for high Mach numbers and provided a data analysis which indicated that such low- speed, low Reynolds 
number corrections could still be applied at high Mach numbers. However, it was necessary to define a 
different Reynolds number. He recommended 8 + , the Reynolds number based on the friction velocity, wall 
density, and boundary-layer thickness. Figure 8, taken from reference 30, shows the domain of importance 
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for including low Reynolds number effects. For values of 8 + below 3000, the effects become more impor- 
tant and, in particular, below 400, they are significant. Lines of constant Mach number indicate that low 
Reynolds number effects can become substantial at high Mach numbers even though the length Reynolds 
number is large. 

Modeling for adverse pressure gradient flows at hypersonic Mach numbers is less well advanced 
because the data base is limited. Figure 9 presents a list of experiments and pertinent test variables. It 
represents a partial, but representative, list of benchmark flows available for model evaluation. Mach 
number is limited to 7 and the wall-to-total-temperature range is mostly adiabatic. For these representative 
flows, eddy viscosity models give reasonably accurate results. Typical comparisons between computations 
and experiments for the skin friction, taken from reference 19, are shown in figure 10. The C-S model with 
the pressure gradient correction (i.e., the p + term) and the higher-order eddy viscosity and Reynolds stress 
models all adequately predict the influence of pressure gradient over a wide range of Reynolds numbers. 


Shock- Wave Boundary-Layer Interaction Flows 

In this section we discuss several examples of shock- wave boundary-layer interaction flows, some 
of which are separated. Figure 11, taken from reference 21, summarizes the status of experiments and 
computation for a variety of compressible flows. We will discuss a limited number of these flows consisting 
of both supersonic and hypersonic cases. 

Figure 12 illustrates the experimental geometry of two hypersonic flows to be discussed. The first is 
a Mach 7 flow about an ogive-cylinder geometry. Two subcases of this flow will be considered: the first is 
the flow over the clean body from the nose rearward, and the second consists of the shock-wave boundary- 
layer interaction on the cylinder produced by a 15° ring-shock generator. Figure 13 shows comparisons 
of predictions and measurements of surface properties for the clean-body case (ref. 15). Three predictions 
are shown, one corresponding to laminar flow and the other two corresponding to turbulent flow obtained 
using the C-S and q — w, o models. Transition was enforced in the modeling at a location about 10 cm 
from the nose. It is apparent that both turbulence models accurately predict surface pressure, skin friction, 
and heat transfer (Stanton number) distributions. 

Results of predictions and measurements of the shock-wave boundary-layer interaction flow on the 
cylinder are shown in figure 14. In this case, measurements of surface-pressure, skin-friction, and heat- 
transfer are compared with predictions made with the C-S, B-L, and the three versions of the g — w model. 
It is apparent from these results that both the zero-equation models and the unmodified q — u, a models 
strongly overpredict the peak pressure in the interaction (fig. 14a). This basically is the result of the 
inability of the three models to adequately predict the extent of separation, which is indicated by the plateau 
in the measured pressure distribution ahead of the interaction. Substantial improvement was obtained 
with the q — w model b, which incorporated the compressibility correction. Results obtained with the 
q — u, c model, which incorporated the modification for heat transfer, were similar to those for model b. 

Skin-friction and heat-transfer distributions for this case are shown in figures 14b and 14c. It is ap- 
parent from these results that although the two zero-equation models give reasonable predictions of peak 
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heating and skin-friction, their predictions in the region of separation are less accurate. Predictions made 
with the unmodified q — u, a model indicate no separation at all and, as a result, grossly overpredict both 
peak heating and skin-friction. The computation made with the q — w , 6 model shows an improvement in 
skin-friction prediction, but still strongly overpredicts peak heating. Finally, the prediction made with the 
q — u>,c model produced results that were in reasonably good agreement with the measurements. 

Computed and measured pressure contours of the flow are shown in figures 15 and 16. Figure 15 
shows measured pressure contours compared with contours obtained with the q — w, c model. It is appar- 
ent that the overall features of the flow are well predicted by the model. Figure 16 compares predictions 
obtained with the q — u, c and the B-L models. This result illustrates that the differences in model predic- 
tions of surface characteristics are accurately reflected in predictions of flow-field variables as well. 

Calculations and measurements of the 7 .5° shock generator case are also discussed in reference 15. 
In this case the flow is attached and the predictions of the two zero-equation models and the q — u, c model 
give similar results that are in good agreement with the experiment. 

The second flow is a compression comer flow, also illustrated in figure 12. Results of measurements 
and computation are shown in figures 17 and 18, corresponding to attached and separated cases, respec- 
tively (ref. 15). Calculations of surface pressure and heat transfer for the attached flow case (fig. 17), made 
with the C-S, B-L, and q — u,c models, indicate reasonably good agreement between computation and ex- 
perimentation. In the separated case (fig. 18) the predictions are also in reasonable agreement, although 
the pressure plateau and extent of separation predicted by the q — u model is better than that predicted 
by the B-L model. In the reattachment zone, both models underpredict overshoots in measured pressure 
and heat-transfer distributions. In this case, computations with the C-S model were unreliable because of 
difficulties in computing boundary-layer and displacement thickness distributions, and are therefore not 
shown. 

It should be noted that the modifications made to the q — w model were general in the sense that no 
arbitrary constants were introduced and then adjusted to improve predictions. Furthermore, the modifica- 
tions introduced did not interfere with or change predictions of simple attached flows (e.g., the clean-body 
flow). This is the type of modification one seeks when improving turbulence models for complex flows. 

The final shock- wave boundary-layer interaction to be discussed is a Mach-3 compression comer flow 
illustrated in figure 19. Calculations of this flow with a comer angle of 20 ° are compared with results of the 
experiment in figure 20, and are discussed in greater detail in reference 7. The turbulence model used in this 
case was the Jones-Launder k — e model. Two wall treatments were investigated, namely, the integration-to- 
the-wall and wall-function procedures. In this (20 °) case, the flow was mildly separated. When predictions 
and measurements are compared, it is clear that noticeable differences in predictions result from different 
wall treatments using the same model. From both skin-friction and pressure distributions it is clear that the 
wall function treatment gives better predictions of separation and surface pressure. In addition, it also gives 
much better agreement with downstream skin-friction distributions than does the integration-to-the-wall 
procedure. Results similar to these were also observed in the 16° and 24° cases (ref. 7). 

The primary reason for the differences between model predictions in this case lies in the low Reynolds 
number (damping) terms of the Jones-Launder model that strongly influence results when the integration- 
to-the-wall procedure is used, but that are inactive when the wall function procedure is used. Although the 
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low Reynolds number terms can produce accurate results for zero pressure gradient attached flows, it is 
apparent that unless they are chosen carefully they can lead to unreliable predictions of complex flows. 


Compressible Shear Layers 

Knowledge of the physics of high-speed shear layers is limited at present. Experiments have shown 
that the far-field spreading angle, a measure of mixing, is reduced considerably, compared to that for 
incompressible flows. Figure 21 illustrates the status. The inverse of the spreading angles obtained from 
various experiments on single- stream mixing layers are shown as a function of Mach number. Reduction 
in spreading angle by a factor of 3 occurs at Mach 5. Various postulates to explain this reduced mixing 
have been proposed, but experimental evidence to substantiate them is lacking. 

Application of incompressible turbulence models, extended to account for compressibility as de- 
scribed in previous sections, fails to produce accurate predictions of the spreading rate. The line labeled 
k — e represents such a prediction. In reference 31, Dash et al. proposed a new model to account for 
compressibility in free shear layers and predicted the spreading angle as shown by the symbols labeled 
k — e, cc. The compressibility correction (cc) involved the empirical function K{M T ) , shown in figure 22, 
where M r is the ratio of the square root of the turbulent kinetic energy to the speed of sound, k x/2 /a, at 
the point of maximum k in the layer. * 

It is noteworthy that the modification also gives reasonable predictions for two-stream supersonic 
mixing. Figure 23 shows a comparison of k — e model predictions compared with the measurements of 
Chinzie (ref. 32). The inverse of the two-stream spreading rate, scaled by the spreading rate <t 0 for which 
one stream is stationary, compares reasonably well with the k — e, cc model prediction. Deviation of the 
data from the modified model prediction for the higher second-stream Mach numbers may be a result of 
free-stream turbulence present in the experiment. 

Vandromme (ref. 17) also reported successful predictions of the single-stream spreading rate. As 
mentioned in the section on modeling, he used the ideas of Rubesin to make compressibility corrections 
to the turbulent kinetic energy and dissipation equations of the k — e model. The results of his predictions 
are shown compared with results of experiments in figure 24. Substantial agreement was achieved. 

More work will be necessary before the compressible mixing layer problem can be considered solved. 
Current modeling modifications are, to a considerable extent, ad hoc and have not been verified for a wide 
range of cases. Futhermore, they are not based on an understanding of the physical mechanisms involved. 
To understand these mechanisms, more experimentation is needed. It should also be noted that research 
is underway at Ames Research Center to use full simulations of compressible shear layers using the time- 
dependent Navier-Stokes equations to provide more complete information on mixing phenomena. It is 
hoped that this research will lead to improved modeling of compressible shear flows. 
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CONCLUDING REMARKS 


In the preceding paragraphs we have described the development and status of turbulence models used 
in the numerical simulation of complex hypersonic flows. In our discussion we emphasized eddy viscosity 
models which constitute the simplest but most widely used class of turbulence models. Two subgroups 
of models were discussed- zero-equation and two-equation models. Each of these models has theoretical 
advantages over the other. For example, two-equation models provide a more general specification of 
turbulent length and velocity scales than zero-equation models, but they often display numerical stability 
problems which are not common to zero-equation models. 

The basic models discussed are similar to those developed originally for incompressible flow. This 
is because in many applications, especially to simple attached boundary-layer flows at low to moderate 
supersonic speeds, the incompressible forms give satisfactory results. As discussed in the text, however, 
there is evidence that these incompressible forms become unsatisfactory as the flow complexity and/or the 
Mach number increase. With respect to flow complexity, it was shown that compressibility corrections 
were necessary to give satisfactory predictions of several hypersonic shock-wave boundary-layer interac- 
tion flows. With respect to Mach number, it was shown that incompressible model forms are unsatisfactory 
for compressible free-shear flows. In this case, too, compressibility corrections could be found which lead 
to satisfactory predictions. 

The status of turbulence modeling for hypersonic flows is still far from complete, however. More 
experimental data and computational comparisons will be necessary to verify and establish the compress- 
ibility corrections made to date. In addition, more experimental and computational work will be needed, 
especially at low Reynolds numbers, because this flow regime is more prevalent at hypersonic speeds, and 
because the available data base in this case is still quite limited. 
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Fig. 1. Composite sketch of a turbulent boundary layer. 



Fig. 2. A log-law representation of hypersonic boundary layer profiles. 



MAISE & MCDONALD M= 5 



Fig. 3. A comparison of the best estimate of shear stress data with incompressible 
measurements and with Maise and McDonald’s compressible approximation. 



22 



DISTANCE FROM WALL, y/5 

(a) Mixing length. 

Fig. 4. Mixing length and scaled eddy viscosity from Maise and McDonald. 
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(b) Scaled eddy viscosity. 
Fig. 4. Concluded. 
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FLAT PLATE SKIN FRICTION 


VAN DRIEST 

COMPUTATIONS BY RUBESIN ET AL. 




Fig. 5. Ability of turbulence models to predict compressibility effects: (a) adiabatic 
wall temperature; (b) Mach 5 and variable wall temperature. 
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Fig. 7. Error in skin friction prediction for two ranges of Reynolds numbers using 
a turbulence model with and without correction for density fluctuations. 
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Fig. 8. Increasing importance of Low Reynolds number effects with Mach number. 
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(a) Transonic flows with impinging shock waves. 
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(b) Supersonic flows with impinging shock waves. 
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(c) Supersonic flows with compression surfaces. 




00 

O) 


<0 

S 

LU 

z 

* 

LU 

u 

Z 

LU 

oc 


o 

a 

a 

oc 

O 

UL 

z 

< 


CO 
-J 
LU 

D 

o 
2 
Z 
O 

P “ 

= I 
2 £ 

6 CC 
> CO 


a 

z 

< 


< 

CC 

o 


Z a 

o - 

6 

CC -I 
LU LU 
Nl 00 
(O -O 


co 

l 

O 

o 

X! 

oo 

W) 

g 

'3 

c 

u 


co 

> 

O 

C 


a 

o 


<D 

a 

D 

C/5 



£ 

60 


Cu 


% 


% 


34 
















15 


10 



Q. 


5 


0 

10 


a 5 


o 

-4 


Fig. 17. Compr 
surface heat transfer 






40 





EXPERIMENT: SETTLES et al. 


Moo = 2.85 Reg ~ 1.6 X 10 6 

TJTj ~ 1 5 0 ~ 2.5 cm 




Fig. 19. Geometry and conditions of a compression corner experiment. 
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TURBULENCE MODEL COMPARISON WITH WALL FUNCTIONS 
AND WITH INTEGRATION TO WALL 


20° RAMP ANGLE 


SURFACE PRESSURE 





SKIN FRICTION 


M^ = 2.79 

O EXPERIMENT, SETTLES etal. 

COMPUTATIONS, k - e MODEL 

WALL FUNCTIONS, yX * 55 

INTEGRATION TO WALL 


Fig. 20. Comparison of computations using k - e model with experiment. 








MACH NUMBER, M 1 

Fig. 21. Spreading rates for compressible 2-D shear layers. Data fits from NASA 
SP-321, 1972. 
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Fig. 22. Compressibility correction for k - e model proposed by Dash, et. al 31 . 








IWNSA 

Natoul A**oi Mutes and 
S|\h.>* Admtmslr.rtion 


1. Report No. 

NASA TM-101079 


4. Title and Subtitle 


Turbulence Modeling for Hypersonic Flows 


Report Documentation Page 


2. Government Accession No. 


3. Recipient's Catalog No. 


5. Report Date 

June 1989 


6. Performing Organization Code 


7. Author(s) 

J. G. Marvin and T. J. Coakley 


9. Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


8. Performing Organization Report No. 

A-89060 

10. Work Unit No. 

505-60-1 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 


15. Supplementary Notes 

Point of Contact: Joe Marvin, Ames Research Center, MS 229-1 , Moffett Field, CA 94035 
(415) 694-5390 or FTS 464-5390 


16. Abstract 

Turbulence modeling for high speed compressible flows is described and discussed. Starting with 
the compressible Navier-Stokes equations, methods of statistical averaging are described by means of 
which the Reynolds-averaged Navier-Stokes equations are developed. Unknown averages in these equa- 
tions are approximated using various closure concepts. Zero-, one-, and two-equation eddy viscosity 
models, algebraic stress models and Reynolds stress transport models are discussed. Computations of 
supersonic and hypersonic flows obtained using several of the models are discussed and compared with 
experimental results. Specific examples include attached boundary layer flows, shock wave boundary 
layer interactions and compressible shear layers. From these examples, conclusions regarding the status 
of modeling and recommendations for future studies are discussed. 


17. Key Words (Suggested by Author(s)) 
Hypersonic 
Boundary layer 
Turbulence modeling 
Shock, Boundary-Layer interaction 


18. Distribution Statement 

Unclassified — Unlimited 


Subject Category — 34 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of pages 

22. Price 

Unclassified 

Unclassified 

47 

A03 


NASA FORM 1626 OCT 86 


For sale by the National Technical Information Service Springfield Virginia 22161 











